Baby Mental Life: Study 2

Exploratory analyses

“Baby Mental Life: Study 2” was conducted on MTurk on 2018-08-04.

Our planned sample was 300 participants, and we anticipated that roughly 80% of recruited participants would pass all of our attention checks, so we initially recruited 378 participants (on the idea that ~80% of 378 ~ 300 participants; note that for administrative purposes we need to recuit participants in batches that were divisible by 9). After filtering out participants who failed at least one of our attention checks, we ended up retaining fewer than 300 participants, so we recruited an additional 16 participants for a total of 394 people recruited. At each stage, we recruited women and men through separate studies, in hopes of acquiring a roughly equal split between genders.

In the end, we ended up with a sample of 304 participants who passed our attention checks, 237 of whom came from unique GPS coordinates.

For this first pass, these data exclude participants where there is another participant with an identical set of GPS coordinates as recorded by Qualtrics.

Each participant assessed children’s mental capacities at 13 target ages between the ages of 0 and 5 years. For each target, they rated 20 mental capacities on a scale from 0 (not at all capable) to 100 (completely capable).

For more details about the study, see our preregistration here.

Here we run some exploratory analyses on these data.

# load required libraries
library(tidyverse)
library(langcog) # source: https://github.com/langcog/langcog-package
library(psych)
library(lme4)
# set theme for ggplots
theme_set(theme_bw())

# run source code (extra home-made functions)
source("./scripts/max_factors_efa.R")
source("./scripts/plot_fun.R")
source("./scripts/reten_fun.R")
source("./scripts/table_fun.R")
source("./scripts/data_prep.R")

NAs introduced by coercionattributes are not identical across measure variables;
they will be droppedJoining, by = "question_qualtrics"

Treating factors as categories

Preliminaries

First, I’ll do the factor analysis, and get a list of the top 5 items by factor:

# load in S1 efa in case we need it
efa_S1 <- readRDS("../study 1/s1_efa.rds")
demo_S1 <- read.csv("../study 1/s1_demo.csv")

# conduct S2 efa
efa_S2 <- fa(d_all, nfactors = 4, rotate = "oblimin", fm = "minres",
             scores = "tenBerge", impute = "median")

Loading required namespace: GPArotation

table_fun(efa_S2, num_items = 5, pos_abs = "abs")

capacity	loading
Factor 1
feeling_distressed	0.82
feeling_overwhelmed	0.82
feeling_frustrated	0.78
feeling_helpless	0.76
feeling_lonely	0.60
Factor 2
having_self_control	0.96
controlling_their_emotions	0.93
telling_right_from_wrong	0.93
planning	0.89
reasoning_about_things	0.89
Factor 3
getting_hungry	0.83
feeling_pain	0.79
feeling_tired	0.67
hearing_sounds	0.58
feeling_physically_uncomfortable	0.49
Factor 4
feeling_happy	0.82
finding_something_funny	0.81
feeling_excited	0.79
loving_somebody	0.64
learning_from_other_people	0.51

Now I’ll use this to define categories of mental capacities:

# get items by factor
factors_S2 <- efa_S2$loadings[] %>%
  data.frame() %>%
  rownames_to_column("capacity") %>%
  gather(factor, loading, -capacity) %>%
  group_by(capacity) %>%
  top_n(1, loading) %>%
  ungroup() %>%
  select(-loading) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR1" = "Negative emotions",
                                      "MR2" = "Cognition & control",
                                      "MR3" = "Bodily sensations",
                                      "MR4" = "Positive/social emotions"),
         factor = factor(factor))
# make new dataframe
d_cat <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(d_demo %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid)))

Joining, by = "capacity"
Joining, by = "subid"

Check reliability (Cronbach’s alpha):

d_cat %>% 
  filter(factor == "MR1") %>% 
  mutate(subid_target = paste(subid, target, sep = "_")) %>%
  select(subid_target, capacity, response) %>%
  spread(capacity, response) %>%
  column_to_rownames("subid_target") %>%
  cor() %>%
  alpha(title = "MR1: Negative emotions")

Reliability analysis  MR1: Negative emotions  
Call: alpha(x = ., title = "MR1: Negative emotions")

 

 Reliability if an item is dropped:

 Item statistics 

d_cat %>% 
  filter(factor == "MR2") %>% 
  mutate(subid_target = paste(subid, target, sep = "_")) %>%
  select(subid_target, capacity, response) %>%
  spread(capacity, response) %>%
  column_to_rownames("subid_target") %>%
  cor() %>%
  alpha(title = "MR2: Cognition & control")

Reliability analysis  MR2: Cognition & control  
Call: alpha(x = ., title = "MR2: Cognition & control")

 

 Reliability if an item is dropped:

 Item statistics 

d_cat %>% 
  filter(factor == "MR3") %>% 
  mutate(subid_target = paste(subid, target, sep = "_")) %>%
  select(subid_target, capacity, response) %>%
  spread(capacity, response) %>%
  column_to_rownames("subid_target") %>%
  cor() %>%
  alpha(title = "MR3: Bodily sensations")

Reliability analysis  MR3: Bodily sensations  
Call: alpha(x = ., title = "MR3: Bodily sensations")

 

 Reliability if an item is dropped:

 Item statistics 

d_cat %>% 
  filter(factor == "MR4") %>% 
  mutate(subid_target = paste(subid, target, sep = "_")) %>%
  select(subid_target, capacity, response) %>%
  spread(capacity, response) %>%
  column_to_rownames("subid_target") %>%
  cor() %>%
  alpha(title = "MR4: Social abilities & positive emotions")

Reliability analysis  MR4: Social abilities & positive emotions  
Call: alpha(x = ., title = "MR4: Social abilities & positive emotions")

 

 Reliability if an item is dropped:

 Item statistics 

And get an average “score” for each of these factors for each participant:

d_cat_scored <- d_cat %>%
  group_by(subid, parent, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

Plots

# bootstrapped means and 95% CIs
d_cat_scored_boot <- d_cat_scored %>%
  group_by(target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_scored,
       aes(x = target_num, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)")

ggplot(d_cat_scored,
       aes(x = target_num, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root-transformed",
       x = "age after square-root transformation (months)")

ggplot(d_cat_scored,
       aes(x = target_ord, y = score, color = factor_names)) +
  facet_grid(~ factor_names) +
  geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_scored_boot, color = "black",
            aes(y = mean, group = factor_names)) +
  geom_pointrange(data = d_cat_scored_boot, color = "black", fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set2", guide = "none") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)")

Regressions

# contrasts(d_cat_scored$target_ord) <- contr.poly(13)
# contrasts(d_cat_scored$factor) <- contr.sum(4)
contrasts(d_cat$target_ord) <- contr.poly(13)
contrasts(d_cat$factor) <- contr.sum(4)

r1 <- lmer(response ~ factor * poly(target_num, 3) +
             (1 + factor + target_num | subid),
           d_cat %>%
             mutate(target_num = scale(target_num, center = F)))
summary(r1)

Linear mixed model fit by REML ['lmerMod']
Formula: response ~ factor * poly(target_num, 3) + (1 + factor + target_num |  
    subid)
   Data: d_cat %>% mutate(target_num = scale(target_num, center = F))

REML criterion at convergence: 536724.6

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-5.5360 -0.4298  0.0225  0.5167  5.5499 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 200.49   14.159                          
          factor1     156.13   12.495    0.65                  
          factor2     145.83   12.076   -0.44 -0.59            
          factor3      69.75    8.352   -0.55 -0.34 -0.24      
          target_num   51.83    7.199   -0.64 -0.42  0.50  0.31
 Residual             334.01   18.276                          
Number of obs: 61620, groups:  subid, 237

Fixed effects:
                               Estimate Std. Error t value
(Intercept)                     68.3739     0.7634  89.569
factor1                          6.9512     0.8216   8.461
factor2                        -38.7255     0.7947 -48.728
factor3                         25.1428     0.5573  45.116
poly(target_num, 3)1          3241.2664    88.4095  36.662
poly(target_num, 3)2         -1237.5812    18.2759 -67.717
poly(target_num, 3)3           532.0777    18.2759  29.114
factor1:poly(target_num, 3)1  -544.0136    31.6548 -17.186
factor2:poly(target_num, 3)1  2670.8981    31.6548  84.376
factor3:poly(target_num, 3)1 -2550.8843    31.6548 -80.584
factor1:poly(target_num, 3)2    32.5828    31.6548   1.029
factor2:poly(target_num, 3)2    38.2604    31.6548   1.209
factor3:poly(target_num, 3)2   881.6516    31.6548  27.852
factor1:poly(target_num, 3)3    48.0128    31.6548   1.517
factor2:poly(target_num, 3)3  -391.9779    31.6548 -12.383
factor3:poly(target_num, 3)3  -354.1489    31.6548 -11.188

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

This is weirdly hard to model… lots of alternative models didn’t converge. Also, those extremely high t-values make me a little suspicious, as do the identical standard errors for the last 9 terms. So, let’s take this with a bit of a grain of salt for now.

Parents vs. non-parents

Study 2

Preliminaries

d_parent <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(d_demo %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid)))

Joining, by = "subid"

parent_counts <- d_demo %>% 
  distinct(ResponseId, Parent) %>% 
  rename(subid = ResponseId, parent = Parent) %>%
  count(parent) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(parent_counts, digits = 3)

parent	n	proportion
No	146	0.616
Yes	88	0.371
NA	3	0.013

d_demo %>%
  count(Parent, ChildrenOldestAge_collapse)

Given these numbers, I think our best bet it just to look at parents vs. non-parents, and not try to separate out parents of young children (too few!).

scores_S2 <- efa_S2$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(d_parent %>% distinct(subid, parent)) %>%
  left_join(d_cat %>% distinct(factor, factor_names)) %>%
  mutate(factor = factor(factor))

Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector

Plots

Factor scores

# bootstrapped means and 95% CIs
d_parent_boot <- scores_S2 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5, # note: too close to dodge
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Exact age, untransformed",
       color = "parent status: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent),
            position = position_dodge(width = 0.3)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.3)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Exact age, square-root transformed",
       color = "parent status: ",
       x = "age after square-root transformation (months)", 
       y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_ord, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot, aes(y = mean, group = parent),
            position = position_dodge(width = 0.5)) +
  geom_pointrange(data = d_parent_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.5)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status",
       subtitle = "Age as ordinal variable",
       color = "parent status: ",
       x = "age (ordinal)", 
       y = "score (NOTE: scales differ)")

Category ‘summary scores’

# bootstrapped means and 95% CIs
d_cat_parent_scored_boot <- d_cat_scored %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_parent_scored_boot,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_parent_scored_boot,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformation",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_parent_scored_boot,
       aes(x = target_ord, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot, 
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)",
       y = "score (NOTE: scales differ)")

All of these visualizations suggest that, relative to non-parents (n = NA), parents (n = NA) tended to perceive greater abilities, in all domains except for bodily sensations. For both negative emotions and positive/social emotions, this difference is greatest in the mid-range of the target ages (12-24 months); for cognition & control, it’s greater at the older end of the target ages (e.g., 3-5 years). It doesn’t look like a huge effect, but it’s intruiging.

Regressions

There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).

contrasts(scores_S2$parent) <- contr.treatment(2, base = 1) # baseline: NON-parents
contrasts(scores_S2$factor) <- contr.sum(4)

r2 <- lmer(score ~ target_num * factor * parent
           + (target_num + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r2)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * parent + (target_num + factor |  
    subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 22480.5

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-9.3553 -0.4388  0.0524  0.5114  5.2974 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 0.3491   0.5908                          
          target_num  0.0165   0.1284   -0.73                  
          factor1     0.1947   0.4413    0.49 -0.40            
          factor2     0.2303   0.4799   -0.66  0.56 -0.40      
          factor3     0.2605   0.5104    0.24  0.00 -0.22 -0.61
 Residual             0.2919   0.5403                          
Number of obs: 12168, groups:  subid, 234

Fixed effects:
                            Estimate Std. Error t value
(Intercept)                -0.424399   0.049602  -8.556
target_num                  0.279097   0.011334  24.624
factor1                     0.078297   0.039260   1.994
factor2                    -0.273518   0.042253  -6.473
factor3                     0.253010   0.044634   5.669
parent2                     0.070601   0.080884   0.873
target_num:factor1         -0.058596   0.006811  -8.603
target_num:factor2          0.190717   0.006811  28.001
target_num:factor3         -0.170985   0.006811 -25.104
target_num:parent2          0.003333   0.018482   0.180
factor1:parent2             0.005619   0.064020   0.088
factor2:parent2            -0.051089   0.068900  -0.741
factor3:parent2             0.031692   0.072783   0.435
target_num:factor1:parent2 -0.002939   0.011107  -0.265
target_num:factor2:parent2  0.053753   0.011107   4.840
target_num:factor3:parent2 -0.046204   0.011107  -4.160

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

The last three coefficients here are where the action is: According to this model, the difference between parents and non-parents in perceived developmental trajectories (target_num:parent2, which is not significant collapsing across factors) is exaggerated in the domain of cognition and control (target_num:factor2:parent2) and diminished in the domain of bodily sensations (target_num:factor3:parent2).

Let’s try adding polynomial effects:

r3 <- lmer(score ~ poly(target_num, 3) * factor * parent
           + (poly(target_num, 1) + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r3)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * parent + (poly(target_num,  
    1) + factor | subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 19808.5

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-10.3042  -0.4652   0.0017   0.5305   5.8317 

Random effects:
 Groups   Name                Variance Std.Dev. Corr                   
 subid    (Intercept)           0.2278  0.4773                         
          poly(target_num, 1) 519.8012 22.7991  -0.51                  
          factor1               0.1982  0.4452   0.45 -0.39            
          factor2               0.2338  0.4835  -0.60  0.54 -0.40      
          factor3               0.2640  0.5138   0.30  0.00 -0.22 -0.61
 Residual                       0.2316  0.4812                         
Number of obs: 12168, groups:  subid, 234

Fixed effects:
                                       Estimate Std. Error t value
(Intercept)                           -0.030561   0.039884  -0.766
poly(target_num, 3)1                  48.855288   1.983984  24.625
poly(target_num, 3)2                 -19.550418   0.613116 -31.887
poly(target_num, 3)3                   8.448995   0.613116  13.780
factor1                               -0.004388   0.038065  -0.115
factor2                               -0.004395   0.041145  -0.107
factor3                                0.011731   0.043587   0.269
parent2                                0.075304   0.065038   1.158
poly(target_num, 3)1:factor1         -10.257028   1.061949  -9.659
poly(target_num, 3)2:factor1           3.413024   1.061949   3.214
poly(target_num, 3)3:factor1          -1.785482   1.061949  -1.681
poly(target_num, 3)1:factor2          33.384462   1.061949  31.437
poly(target_num, 3)2:factor2           3.384910   1.061949   3.187
poly(target_num, 3)3:factor2          -6.433547   1.061949  -6.058
poly(target_num, 3)1:factor3         -29.930494   1.061949 -28.184
poly(target_num, 3)2:factor3           8.594647   1.061949   8.093
poly(target_num, 3)3:factor3          -3.251241   1.061949  -3.062
poly(target_num, 3)1:parent2           0.583455   3.235228   0.180
poly(target_num, 3)2:parent2          -2.031708   0.999792  -2.032
poly(target_num, 3)3:parent2           1.883539   0.999792   1.884
factor1:parent2                        0.001472   0.062072   0.024
factor2:parent2                        0.024763   0.067094   0.369
factor3:parent2                       -0.033507   0.071076  -0.471
poly(target_num, 3)1:factor1:parent2  -0.514473   1.731690  -0.297
poly(target_num, 3)2:factor1:parent2  -0.041808   1.731690  -0.024
poly(target_num, 3)3:factor1:parent2   1.099314   1.731690   0.635
poly(target_num, 3)1:factor2:parent2   9.409337   1.731690   5.434
poly(target_num, 3)2:factor2:parent2  -2.361622   1.731690  -1.364
poly(target_num, 3)3:factor2:parent2  -1.057673   1.731690  -0.611
poly(target_num, 3)1:factor3:parent2  -8.087893   1.731690  -4.671
poly(target_num, 3)2:factor3:parent2   6.669873   1.731690   3.852
poly(target_num, 3)3:factor3:parent2  -4.830063   1.731690  -2.789

Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

A lot to sort through here, but here are some observations:

• Collapsing across factors, it seems like parents might perceive exaggerated non-linearities? (poly(target_num, 3)2:parent2 and poly(target_num, 3)3:parent2)
• Compared to collapsing across factors, not a lot of additional parent vs. non-parent differences in negative emotions (factor 1)
• Compared to collapsing across factors, parents perceived more dramatic (linear) growth in cognition/control (factor 2; poly(target_num, 3)1:factor2:parent2)
• Compared to collapsing across factors, in the domain of bodily sensations (factor 3), parents perceived less growth (poly(target_num, 3)1:factor3:parent2). I’m not sure how to interpret the differences in non-linearity here.

You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.

Study 1

Preliminaries

parent_counts_S1 <- demo_S1 %>% 
  distinct(ResponseId, Parent) %>% 
  rename(subid = ResponseId, parent = Parent) %>%
  count(parent) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(parent_counts_S1, digits = 3)

parent	n	proportion
No	147	0.653
Yes	77	0.342
NA	1	0.004

scores_S1 <- efa_S1$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(demo_S1 %>% distinct(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid))) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR4" = "Negative emotions (S1 F4)",
                                      "MR1" = "Cognition & control (S1 F1)",
                                      "MR2" = "Bodily sensations (S1 F2)",
                                      "MR3" = "Positive/social emotions (S1 F3)"),
         factor = factor(factor))

Joining, by = "subid"

Plots

Factor scores

# bootstrapped means and 95% CIs
d_parent_boot_S1 <- scores_S1 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S1 %>%
  filter(!is.na(parent)) %>%
  ggplot(aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_parent_boot_S1,
            position = position_dodge(width = 5),
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_parent_boot_S1, 
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = parent), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by parent status (STUDY 1)",
       subtitle = "Exact age, untransformed",
       color = "parent status: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

I won’t bother to re-plot with different treatments of age, since there are only three target ages here.

Category ‘summary scores’

# make new dataframe
d_cat_S1 <- d_all_S1 %>%
  rename(subid_target = X) %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(demo_S1 %>%
              select(ResponseId, Parent) %>%
              rename(subid = ResponseId, parent = Parent) %>%
              mutate(subid = as.character(subid))) %>%
  filter(!is.na(factor)) # drop extra items from S1 not in S2

Joining, by = "capacity"
Joining, by = "subid"

And get an average “score” for each of these factors for each participant:

d_cat_scored_S1 <- d_cat_S1 %>%
  group_by(subid, parent, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

# bootstrapped means and 95% CIs
d_cat_parent_scored_boot_S1 <- d_cat_scored_S1 %>%
  filter(!is.na(parent)) %>%
  group_by(parent, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_parent_scored_boot_S1,
       aes(x = target_num, y = score, color = parent)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_parent_scored_boot_S1,
            position = position_dodge(width = 5),
            aes(y = mean, group = parent)) +
  geom_pointrange(data = d_cat_parent_scored_boot_S1,
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

Looking quickly at these two visualizations of Study 1, I’d say that there’s a hint of the same patterns in the domains of negative emotions and positive social emotions (parents attributed more than non-parents, especially at 9 months), and in cognition/control (parents attributed more than non-parents at 5 years). So, generally consistent?

I won’t run regression analyses just now.

Men vs. women

Study 2

Preliminaries

d_gender <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(d_demo %>%
              select(ResponseId, GenderSex) %>%
              rename(subid = ResponseId, gender = GenderSex) %>%
              mutate(subid = as.character(subid)))

Joining, by = "subid"

gender_counts <- d_demo %>% 
  distinct(ResponseId, GenderSex) %>% 
  rename(subid = ResponseId, gender = GenderSex) %>%
  count(gender) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(gender_counts, digits = 3)

gender	n	proportion
Female	123	0.519
Male	114	0.481

scores_S2 <- efa_S2$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(d_gender %>% distinct(subid, gender)) %>%
  left_join(d_cat %>% distinct(factor, factor_names)) %>%
  mutate(factor = factor(factor), gender = factor(gender))

Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector

Plots

Factor scores

# bootstrapped means and 95% CIs
d_gender_boot <- scores_S2 %>%
  filter(!is.na(gender)) %>%
  group_by(gender, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S2 %>%
  filter(!is.na(gender)) %>%
  ggplot(aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_gender_boot, aes(y = mean, group = gender)) +
  geom_pointrange(data = d_gender_boot, fatten = 1.5, # note: too close to dodge
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = gender), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant gender",
       subtitle = "Exact age, untransformed",
       color = "gender: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(gender)) %>%
  ggplot(aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_gender_boot, aes(y = mean, group = gender),
            position = position_dodge(width = 0.3)) +
  geom_pointrange(data = d_gender_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.3)) +
  # geom_smooth(aes(group = gender), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant gender",
       subtitle = "Exact age, square-root transformed",
       color = "participant gender: ",
       x = "age after square-root transformation (months)", 
       y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(gender)) %>%
  ggplot(aes(x = target_ord, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_gender_boot, aes(y = mean, group = gender),
            position = position_dodge(width = 0.5)) +
  geom_pointrange(data = d_gender_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.5)) +
  # geom_smooth(aes(group = gender), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant gender",
       subtitle = "Age as ordinal variable",
       color = "participant gender: ",
       x = "age (ordinal)", 
       y = "score (NOTE: scales differ)")

Category ‘summary scores’

# bootstrapped means and 95% CIs
d_cat_gender_scored_boot <- d_cat_scored %>%
  left_join(d_demo %>% distinct(ResponseId, GenderSex) %>%
              mutate(ResponseId = as.character(ResponseId)) %>%
              rename(subid = ResponseId, gender = GenderSex)) %>%
  filter(!is.na(gender)) %>%
  group_by(gender, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

Joining, by = "subid"

ggplot(d_cat_gender_scored_boot,
       aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot, 
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_gender_scored_boot,
       aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot, 
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformation",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_gender_scored_boot,
       aes(x = target_ord, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot, 
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)",
       y = "score (NOTE: scales differ)")

All of these visualizations suggest that, relative to women, men tended to perceive fewer abilities in the negative emotions domain (across the age range) and in the bodily sensations domain (especially for very young infants), and perhaps in the positive/social emotions domain (for toddlers). In contrast, if anythign men seemed to perceive more cognition & control abilities in the middle age ranges (~ 1 year). None of these are huge effects, but they are intriguing.

Regressions

There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).

contrasts(scores_S2$gender) <- contr.treatment(2, base = 1) # baseline: FEMALE
contrasts(scores_S2$factor) <- contr.sum(4)

r2 <- lmer(score ~ target_num * factor * gender
           + (target_num + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r2)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * gender + (target_num + factor |  
    subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 22684.2

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-9.3855 -0.4361  0.0520  0.5098  5.3101 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 0.34624  0.5884                          
          target_num  0.01648  0.1284   -0.73                  
          factor1     0.18436  0.4294    0.49 -0.41            
          factor2     0.22582  0.4752   -0.65  0.57 -0.39      
          factor3     0.25876  0.5087    0.24  0.00 -0.23 -0.61
 Residual             0.29011  0.5386                          
Number of obs: 12324, groups:  subid, 237

Fixed effects:
                            Estimate Std. Error t value
(Intercept)                -0.352499   0.053820  -6.550
target_num                  0.279569   0.012339  22.658
factor1                     0.198080   0.041761   4.743
factor2                    -0.398918   0.045618  -8.745
factor3                     0.304166   0.048464   6.276
gender2                    -0.085230   0.077601  -1.098
target_num:factor1         -0.073166   0.007398  -9.890
target_num:factor2          0.232390   0.007398  31.413
target_num:factor3         -0.205753   0.007398 -27.813
target_num:gender2         -0.001484   0.017791  -0.083
factor1:gender2            -0.236684   0.060213  -3.931
factor2:gender2             0.211000   0.065775   3.208
factor3:gender2            -0.080956   0.069878  -1.159
target_num:factor1:gender2  0.028011   0.010667   2.626
target_num:factor2:gender2 -0.044941   0.010667  -4.213
target_num:factor3:gender2  0.037002   0.010667   3.469

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

The last six coefficients here are where the action is: According to this model, the difference between male and female in overall mental capacity attributions (gender2, which is not significant collapsing across factors and targets) is exaggerated in the domain of negative emotions (factor1:gender2), diminished in the domain of cognition and control (factor2:gender2), and neither diminshed nor exaggerated in the domain of bodily sensations (factor3:gender2). Meanwhile, the difference between male and female in perceived developmental trajectories (target_num:gender2, which is not significant collapsing across factors) is (likewise) exaggerated in the domain of negative emotions (target_num:factor1:gender2), diminished in the domain of cognition and control (target_num:factor2:gender2) and exaggerated in the domain of bodily sensations (target_num:factor3:gender2).

Let’s try adding polynomial effects:

r3 <- lmer(score ~ poly(target_num, 3) * factor * gender
           + (poly(target_num, 1) + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r3)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * gender + (poly(target_num,  
    1) + factor | subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 20032.2

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-10.2922  -0.4698   0.0025   0.5288   5.8429 

Random effects:
 Groups   Name                Variance Std.Dev. Corr                   
 subid    (Intercept)           0.2251  0.4745                         
          poly(target_num, 1) 519.0181 22.7820  -0.51                  
          factor1               0.1878  0.4333   0.45 -0.40            
          factor2               0.2292  0.4788  -0.59  0.56 -0.39      
          factor3               0.2622  0.5120   0.29  0.00 -0.23 -0.61
 Residual                       0.2312  0.4809                         
Number of obs: 12324, groups:  subid, 237

Fixed effects:
                                      Estimate Std. Error t value
(Intercept)                            0.04200    0.04320   0.972
poly(target_num, 3)1                  48.93792    2.15991  22.657
poly(target_num, 3)2                 -19.76380    0.66749 -29.609
poly(target_num, 3)3                   8.89729    0.66749  13.329
factor1                                0.09484    0.04043   2.345
factor2                               -0.07099    0.04441  -1.599
factor3                                0.01383    0.04733   0.292
gender2                               -0.08732    0.06229  -1.402
poly(target_num, 3)1:factor1         -12.80752    1.15613 -11.078
poly(target_num, 3)2:factor1           2.82583    1.15613   2.444
poly(target_num, 3)3:factor1          -0.97761    1.15613  -0.846
poly(target_num, 3)1:factor2          40.67937    1.15613  35.186
poly(target_num, 3)2:factor2           3.21977    1.15613   2.785
poly(target_num, 3)3:factor2          -7.38714    1.15613  -6.390
poly(target_num, 3)1:factor3         -36.01663    1.15613 -31.153
poly(target_num, 3)2:factor3          13.88965    1.15613  12.014
poly(target_num, 3)3:factor3          -7.53568    1.15613  -6.518
poly(target_num, 3)1:gender2          -0.25972    3.11428  -0.083
poly(target_num, 3)2:gender2          -0.60651    0.96243  -0.630
poly(target_num, 3)3:gender2           0.27380    0.96243   0.284
factor1:gender2                       -0.19716    0.05830  -3.382
factor2:gender2                        0.14758    0.06403   2.305
factor3:gender2                       -0.02874    0.06824  -0.421
poly(target_num, 3)1:factor1:gender2   4.90328    1.66698   2.941
poly(target_num, 3)2:factor1:gender2   0.95398    1.66698   0.572
poly(target_num, 3)3:factor1:gender2  -0.58404    1.66698  -0.350
poly(target_num, 3)1:factor2:gender2  -7.86688    1.66698  -4.719
poly(target_num, 3)2:factor2:gender2  -1.80751    1.66698  -1.084
poly(target_num, 3)3:factor2:gender2   1.33095    1.66698   0.798
poly(target_num, 3)1:factor3:gender2   6.47711    1.66698   3.886
poly(target_num, 3)2:factor3:gender2  -5.73172    1.66698  -3.438
poly(target_num, 3)3:factor3:gender2   5.11585    1.66698   3.069

Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

A lot to sort through here, but here are some observations:

• Most of the differences across genders seem to be in the linear (not quadratic or cubic) components
• The one exception to this is in the bodily sensations domain (factor 3), where men seem to have perceived the shape of the developmental trajectory differently both in terms of its quadratic (poly(target_num, 2)1:factor2:gender2) and its cubic (poly(target_num, 3)1:factor2:gender2) components. But I’m having a hard time interpreting the coefficients at this point.

You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.

Study 1

Preliminaries

gender_counts_S1 <- demo_S1 %>% 
  distinct(ResponseId, GenderSex) %>% 
  rename(subid = ResponseId, gender = GenderSex) %>%
  count(gender) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(gender_counts_S1, digits = 3)

gender	n	proportion
Female	96	0.427
Male	127	0.564
Other (please describe)	1	0.004
Prefer not to say	1	0.004

scores_S1 <- efa_S1$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(demo_S1 %>% distinct(ResponseId, GenderSex) %>%
              rename(subid = ResponseId, gender = GenderSex) %>%
              mutate(subid = as.character(subid))) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR4" = "Negative emotions (S1 F4)",
                                      "MR1" = "Cognition & control (S1 F1)",
                                      "MR2" = "Bodily sensations (S1 F2)",
                                      "MR3" = "Positive/social emotions (S1 F3)"),
         factor = factor(factor))

Joining, by = "subid"

Plots

Factor scores

# bootstrapped means and 95% CIs
d_gender_boot_S1 <- scores_S1 %>%
  filter(!is.na(gender)) %>%
  group_by(gender, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S1 %>%
  filter(gender %in% c("Female", "Male")) %>%
  ggplot(aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_gender_boot_S1 %>%
              filter(gender %in% c("Female", "Male")),
            position = position_dodge(width = 5),
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_gender_boot_S1 %>%
                    filter(gender %in% c("Female", "Male")), 
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = gender), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant gender (STUDY 1)",
       subtitle = "Exact age, untransformed",
       color = "participant gender: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

I won’t bother to re-plot with different treatments of age, since there are only three target ages here.

Category ‘summary scores’

# make new dataframe
d_cat_S1 <- d_all_S1 %>%
  rename(subid_target = X) %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(demo_S1 %>%
              select(ResponseId, GenderSex) %>%
              rename(subid = ResponseId, gender = GenderSex) %>%
              mutate(subid = as.character(subid))) %>%
  filter(!is.na(factor)) # drop extra items from S1 not in S2

Joining, by = "capacity"
Joining, by = "subid"

And get an average “score” for each of these factors for each participant:

d_cat_scored_S1 <- d_cat_S1 %>%
  group_by(subid, gender, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

# bootstrapped means and 95% CIs
d_cat_gender_scored_boot_S1 <- d_cat_scored_S1 %>%
  filter(!is.na(gender)) %>%
  group_by(gender, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_gender_scored_boot_S1 %>%
         filter(gender %in% c("Female", "Male")),
       aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot_S1 %>%
              filter(gender %in% c("Female", "Male")),
            position = position_dodge(width = 5),
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
                    filter(gender %in% c("Female", "Male")),
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_gender_scored_boot_S1 %>%
         filter(gender %in% c("Female", "Male")),
       aes(x = target_num, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot_S1 %>%
              filter(gender %in% c("Female", "Male")),
            position = position_dodge(width = 1),
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
                    filter(gender %in% c("Female", "Male")),
                  position = position_dodge(width = 1),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformed",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_gender_scored_boot_S1 %>%
         filter(gender %in% c("Female", "Male")),
       aes(x = target_ord, y = score, color = gender)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_gender_scored_boot_S1 %>%
              filter(gender %in% c("Female", "Male")),
            position = position_dodge(width = 0.5),
            aes(y = mean, group = gender)) +
  geom_pointrange(data = d_cat_gender_scored_boot_S1 %>%
                    filter(gender %in% c("Female", "Male")),
                  position = position_dodge(width = 0.5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  # scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "age (ordinal)",
       x = "target",
       y = "score (NOTE: scales differ)")

Looking quickly at these two visualizations of Study 1, I’d say that they are similar to the visualizations for Study 2, both in the negative emotions and bodily domains. There a hint of more of a difference here in the positive & social emotions domain. So, generally consistent?

I won’t run regression analyses just now.

Higher vs. lower levels of education

Study 2

Preliminaries

d_edu <- d_all %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(d_demo %>%
              select(ResponseId, Education) %>%
              mutate(edu = case_when(as.numeric(Education) < 8 ~ "less than a BA",
                                     as.numeric(Education) >= 8 ~ "BA or more"),
                     edu = factor(edu, 
                                  levels = c("less than a BA",
                                             "BA or more"))) %>%
              rename(subid = ResponseId) %>%
              mutate(subid = as.character(subid)))

Joining, by = "subid"

edu_counts <- d_demo %>% 
  distinct(ResponseId, Education) %>% 
  rename(subid = ResponseId, edu = Education) %>%
  count(edu) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(edu_counts, digits = 3)

edu	n	proportion
Some high school, no diploma	2	0.008
High school graduate, diploma or equivalent (including GED)	31	0.131
Some college credit, no degree	61	0.257
Trade school, technical school, or vocational school	11	0.046
Associate's degree (for example, AA, AS)	26	0.110
Bachelor's degree (for example, BA, BS)	80	0.338
Master's degree (for example, MA, MS)	25	0.105
Doctor or professional degree (for example, PhD, JD, MD, MBA)	1	0.004

scores_S2 <- efa_S2$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(d_edu %>% distinct(subid, edu)) %>%
  left_join(d_cat %>% distinct(factor, factor_names)) %>%
  mutate(factor = factor(factor), edu = factor(edu))

Joining, by = "subid"
Joining, by = "factor"
Column `factor` joining character vector and factor, coercing into character vector

Plots

Factor scores

# bootstrapped means and 95% CIs
d_edu_boot <- scores_S2 %>%
  filter(!is.na(edu)) %>%
  group_by(edu, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S2 %>%
  filter(!is.na(edu)) %>%
  ggplot(aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_edu_boot, aes(y = mean, group = edu)) +
  geom_pointrange(data = d_edu_boot, fatten = 1.5, # note: too close to dodge
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = edu), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant edu",
       subtitle = "Exact age, untransformed",
       color = "edu: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(edu)) %>%
  ggplot(aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_edu_boot, aes(y = mean, group = edu),
            position = position_dodge(width = 0.3)) +
  geom_pointrange(data = d_edu_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.3)) +
  # geom_smooth(aes(group = edu), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant edu",
       subtitle = "Exact age, square-root transformed",
       color = "participant edu: ",
       x = "age after square-root transformation (months)", 
       y = "score (NOTE: scales differ)")

scores_S2 %>%
  filter(!is.na(edu)) %>%
  ggplot(aes(x = target_ord, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_edu_boot, aes(y = mean, group = edu),
            position = position_dodge(width = 0.5)) +
  geom_pointrange(data = d_edu_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper),
                  position = position_dodge(width = 0.5)) +
  # geom_smooth(aes(group = edu), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant edu",
       subtitle = "Age as ordinal variable",
       color = "participant edu: ",
       x = "age (ordinal)", 
       y = "score (NOTE: scales differ)")

Category ‘summary scores’

# bootstrapped means and 95% CIs
d_cat_edu_scored_boot <- d_cat_scored %>%
  left_join(d_demo %>% distinct(ResponseId, Education) %>%
              mutate(ResponseId = as.character(ResponseId),
                     edu = case_when(as.numeric(Education) < 8 ~ "less than a BA",
                                     as.numeric(Education) >= 8 ~ "BA or more"),
                     edu = factor(edu,
                                  levels = c("less than a BA",
                                             "BA or more"))) %>%
              rename(subid = ResponseId)) %>%
  filter(!is.na(edu)) %>%
  group_by(edu, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

Joining, by = "subid"

ggplot(d_cat_edu_scored_boot,
       aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot, 
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_edu_scored_boot,
       aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot, 
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformation",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_edu_scored_boot,
       aes(x = target_ord, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot, 
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot, fatten = 1.5,
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Age as an ordinal variable",
       x = "age (ordinal)",
       y = "score (NOTE: scales differ)")

All of these visualizations suggest that, relative to less educated participants, more educated participants tended to perceive slightly fewer abilities in the negative emotions domain, but were otherwise generally very similar.

Regressions

There are many ways that we could choose to model this - I’ll try out the factor scores first, adapting our primary regression models (not in this notebook).

contrasts(scores_S2$edu) <- contr.treatment(2, base = 1) # baseline: less educated
contrasts(scores_S2$factor) <- contr.sum(4)

r2 <- lmer(score ~ target_num * factor * edu
           + (target_num + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r2)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ target_num * factor * edu + (target_num + factor | subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 22726.8

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-9.4244 -0.4335  0.0513  0.5090  5.3716 

Random effects:
 Groups   Name        Variance Std.Dev. Corr                   
 subid    (Intercept) 0.34804  0.5899                          
          target_num  0.01646  0.1283   -0.73                  
          factor1     0.19408  0.4406    0.49 -0.40            
          factor2     0.23103  0.4807   -0.66  0.56 -0.41      
          factor3     0.25862  0.5085    0.24  0.00 -0.21 -0.61
 Residual             0.29081  0.5393                          
Number of obs: 12324, groups:  subid, 237

Fixed effects:
                         Estimate Std. Error t value
(Intercept)             -0.391254   0.052284  -7.483
target_num               0.282154   0.011952  23.607
factor1                  0.081491   0.041379   1.969
factor2                 -0.307750   0.044657  -6.891
factor3                  0.290749   0.046956   6.192
edu2                    -0.005014   0.078179  -0.064
target_num:factor1      -0.054535   0.007177  -7.598
target_num:factor2       0.208371   0.007177  29.033
target_num:factor3      -0.194921   0.007177 -27.159
target_num:edu2         -0.007374   0.017871  -0.413
factor1:edu2             0.006130   0.061873   0.099
factor2:edu2             0.023088   0.066774   0.346
factor3:edu2            -0.057068   0.070211  -0.813
target_num:factor1:edu2 -0.011532   0.010732  -1.075
target_num:factor2:edu2  0.005372   0.010732   0.501
target_num:factor3:edu2  0.015575   0.010732   1.451

Correlation matrix not shown by default, as p = 16 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Nada.

Let’s try adding polynomial effects:

r3 <- lmer(score ~ poly(target_num, 3) * factor * edu
           + (poly(target_num, 1) + factor | subid),
           scores_S2 %>%
             mutate(target_num = target_num/12))
summary(r3)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ poly(target_num, 3) * factor * edu + (poly(target_num,  
    1) + factor | subid)
   Data: scores_S2 %>% mutate(target_num = target_num/12)

REML criterion at convergence: 20092.9

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-10.3685  -0.4600   0.0056   0.5311   5.9055 

Random effects:
 Groups   Name                Variance Std.Dev. Corr                   
 subid    (Intercept)           0.2270  0.4764                         
          poly(target_num, 1) 518.4274 22.7690  -0.51                  
          factor1               0.1975  0.4444   0.45 -0.39            
          factor2               0.2344  0.4842  -0.60  0.55 -0.41      
          factor3               0.2620  0.5119   0.29  0.00 -0.22 -0.61
 Residual                       0.2322  0.4818                         
Number of obs: 12324, groups:  subid, 237

Fixed effects:
                                    Estimate Std. Error t value
(Intercept)                         0.006897   0.042031   0.164
poly(target_num, 3)1               49.390345   2.092247  23.606
poly(target_num, 3)2              -21.317213   0.648102 -32.892
poly(target_num, 3)3               10.204882   0.648102  15.746
factor1                             0.004537   0.040120   0.113
factor2                            -0.013717   0.043494  -0.315
factor3                             0.015694   0.045850   0.342
edu2                               -0.015421   0.062848  -0.245
poly(target_num, 3)1:factor1       -9.546136   1.122546  -8.504
poly(target_num, 3)2:factor1        2.830921   1.122546   2.522
poly(target_num, 3)3:factor1       -0.776043   1.122546  -0.691
poly(target_num, 3)1:factor2       36.474744   1.122546  32.493
poly(target_num, 3)2:factor2        3.825457   1.122546   3.408
poly(target_num, 3)3:factor2       -7.834359   1.122546  -6.979
poly(target_num, 3)1:factor3      -34.120438   1.122546 -30.396
poly(target_num, 3)2:factor3       11.638223   1.122546  10.368
poly(target_num, 3)3:factor3       -5.846249   1.122546  -5.208
poly(target_num, 3)1:edu2          -1.290868   3.128486  -0.413
poly(target_num, 3)2:edu2           2.820908   0.969091   2.911
poly(target_num, 3)3:edu2          -2.629115   0.969091  -2.713
factor1:edu2                       -0.010143   0.059991  -0.169
factor2:edu2                        0.030668   0.065035   0.472
factor3:edu2                       -0.035090   0.068559  -0.512
poly(target_num, 3)1:factor1:edu2  -2.018615   1.678516  -1.203
poly(target_num, 3)2:factor1:edu2   1.014589   1.678516   0.604
poly(target_num, 3)3:factor1:edu2  -1.078782   1.678516  -0.643
poly(target_num, 3)1:factor2:edu2   0.940291   1.678516   0.560
poly(target_num, 3)2:factor2:edu2  -3.298146   1.678516  -1.965
poly(target_num, 3)3:factor2:edu2   2.431313   1.678516   1.448
poly(target_num, 3)1:factor3:edu2   2.726340   1.678516   1.624
poly(target_num, 3)2:factor3:edu2  -1.130454   1.678516  -0.673
poly(target_num, 3)3:factor3:edu2   1.724635   1.678516   1.027

Correlation matrix not shown by default, as p = 32 > 12.
Use print(x, correlation=TRUE)  or
     vcov(x)     if you need it

Again, nothing substantial - some indications of some education-related differences in perceptions of non-linear effects? But hard to make too much sense of that.

You could imagine re-running any of these with a square-root transformation on target age, and/or with the ‘summary scores’ by category, but I’m not going to do that right now.

Study 1

Preliminaries

edu_counts_S1 <- demo_S1 %>% 
  distinct(ResponseId, Education) %>% 
  mutate(edu = case_when(
    as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
    as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
         edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
  rename(subid = ResponseId) %>%
  count(edu) %>% 
  mutate(proportion = n/sum(n))
knitr::kable(edu_counts_S1, digits = 3)

edu	n	proportion
less than a BA	132	0.587
BA or more	93	0.413

scores_S1 <- efa_S1$scores %>%
  data.frame() %>%
  rownames_to_column("subid_target") %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(factor, score, -c(subid, starts_with("target"))) %>%
  left_join(demo_S1 %>% distinct(ResponseId, Education) %>%
              mutate(edu = case_when(
                as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
                as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
                edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
              rename(subid = ResponseId) %>%
              mutate(subid = as.character(subid))) %>%
  mutate(factor_names = recode_factor(factor,
                                      "MR4" = "Negative emotions (S1 F4)",
                                      "MR1" = "Cognition & control (S1 F1)",
                                      "MR2" = "Bodily sensations (S1 F2)",
                                      "MR3" = "Positive/social emotions (S1 F3)"),
         factor = factor(factor))

Joining, by = "subid"

Plots

Factor scores

# bootstrapped means and 95% CIs
d_edu_boot_S1 <- scores_S1 %>%
  filter(!is.na(edu)) %>%
  group_by(edu, target, target_num, target_ord, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup() %>%
  data.frame()

scores_S1 %>%
  filter(!is.na(edu)) %>%
  ggplot(aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.05) +
  geom_line(data = d_edu_boot_S1 %>%
              filter(!is.na(edu)),
            position = position_dodge(width = 5),
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_edu_boot_S1 %>%
                    filter(!is.na(edu)), 
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = edu), #color = "black",
  #             method = "lm", formula = "y ~ poly(x, 3)") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(-10, 10, 0.5)) +
  labs(title = "Developmental trajectories of factor scores, by participant education (STUDY 1)",
       subtitle = "Exact age, untransformed",
       color = "participant edu: ",
       x = "target age (months)", y = "score (NOTE: scales differ)")

I won’t bother to re-plot with different treatments of age, since there are only three target ages here.

Category ‘summary scores’

# make new dataframe
d_cat_S1 <- d_all_S1 %>%
  rename(subid_target = X) %>%
  mutate(subid = gsub("_.*$", "", subid_target),
         target = gsub("^.*_", "", subid_target)) %>%
  select(-subid_target) %>%
  mutate(target_num = recode(target,
                             "target00mo" = 0,
                             "target0Xmo" = 4/30,
                             "target01mo" = 1,
                             "target02mo" = 2,
                             "target04mo" = 4,
                             "target06mo" = 6,
                             "target09mo" = 9,
                             "target12mo" = 12,
                             "target18mo" = 18,
                             "target24mo" = 24,
                             "target36mo" = 36,
                             "target48mo" = 48,
                             "target60mo" = 60),
         target_ord = recode_factor(target,
                                    "target00mo" = "newborns",
                                    "target0Xmo" = "4-day-olds",
                                    "target01mo" = "1-month-olds",
                                    "target02mo" = "2-month-olds",
                                    "target04mo" = "4-month-olds",
                                    "target06mo" = "6-month-olds",
                                    "target09mo" = "9-month-olds",
                                    "target12mo" = "12-month-olds",
                                    "target18mo" = "18-month-olds",
                                    "target24mo" = "2-year-olds",
                                    "target36mo" = "3-year-olds",
                                    "target48mo" = "4-year-olds",
                                    "target60mo" = "5-year-olds")) %>%
  gather(capacity, response, -c(subid, starts_with("target"))) %>%
  left_join(factors_S2) %>%
  left_join(demo_S1 %>%
              select(ResponseId, Education) %>%
              mutate(edu = case_when(
                as.numeric(Education) %in% c(1, 4, 6, 7, 8) ~ "less than a BA",
                as.numeric(Education) %in% c(2, 3, 5) ~ "BA or more"),
                edu = factor(edu, levels = c("less than a BA", "BA or more"))) %>%
              rename(subid = ResponseId) %>%
              mutate(subid = as.character(subid))) %>%
  filter(!is.na(factor)) # drop extra items from S1 not in S2

Joining, by = "capacity"
Joining, by = "subid"

And get an average “score” for each of these factors for each participant:

d_cat_scored_S1 <- d_cat_S1 %>%
  group_by(subid, edu, 
           target, target_num, target_ord, 
           factor, factor_names) %>%
  summarise(score = mean(response, na.rm = T)) %>%
  ungroup()

# bootstrapped means and 95% CIs
d_cat_edu_scored_boot_S1 <- d_cat_scored_S1 %>%
  filter(!is.na(edu)) %>%
  group_by(edu, target, target_num, target_ord, factor, factor_names) %>%
  multi_boot_standard("score") %>%
  ungroup()

ggplot(d_cat_edu_scored_boot_S1 %>%
         filter(!is.na(edu)),
       aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot_S1 %>%
              filter(!is.na(edu)),
            position = position_dodge(width = 5),
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
                    filter(!is.na(edu)),
                  position = position_dodge(width = 5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, untransformed",
       x = "target age (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_edu_scored_boot_S1 %>%
         filter(!is.na(edu)),
       aes(x = target_num, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot_S1 %>%
              filter(!is.na(edu)),
            position = position_dodge(width = 1),
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
                    filter(!is.na(edu)),
                  position = position_dodge(width = 1),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Paired") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  scale_x_continuous(breaks = seq(0, 60, 12), trans = "sqrt") +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "Exact age, square-root transformed",
       x = "age after square-root transformation (months)",
       y = "score (NOTE: scales differ)")

ggplot(d_cat_edu_scored_boot_S1 %>%
         filter(edu %in% c("Female", "Male")),
       aes(x = target_ord, y = score, color = edu)) +
  facet_wrap(~ factor_names, ncol = 4, scales = "free") +
  # geom_line(aes(group = subid), alpha = 0.15) +
  geom_line(data = d_cat_edu_scored_boot_S1 %>%
              filter(edu %in% c("Female", "Male")),
            position = position_dodge(width = 0.5),
            aes(y = mean, group = edu)) +
  geom_pointrange(data = d_cat_edu_scored_boot_S1 %>%
                    filter(edu %in% c("Female", "Male")),
                  position = position_dodge(width = 0.5),
                  aes(y = mean, ymin = ci_lower, ymax = ci_upper)) +
  # geom_smooth(aes(group = factor_names),
  #             method = "lm", formula = "y ~ poly(x, 3)",
  #             color = "black") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        legend.position = "bottom") +
  # scale_x_continuous(breaks = seq(0, 60, 12)) +
  scale_y_continuous(breaks = seq(0, 100, 10)) +
  labs(title = "Developmental trajectories of category 'summary scores'",
       subtitle = "age (ordinal)",
       x = "target",
       y = "score (NOTE: scales differ)")

Error: Faceting variables must have at least one value

Looking quickly at these two visualizations of Study 1, I’d say that they are similar to the visualizations for Study 2 - perhaps a little more dramatic, especially in the negative emotions and bodily sensations domains.

I won’t run regression analyses just now.